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A physically  realizable  stationary,  Gaussian, 
random  load  is  simulated  digitally  and  employed  as  the 
forcing  function  in  the  equation  of  motion  of  a damped, 
elastic  beam  whose  resistance  to  deformation  is  due  to 
bending  and  stretching.  Rotations  are  considered  to  be 
finite,  the  squares  of  slopes  being  negligible  compared 
with  unity,  but  of  the  order  of  magnitude  of  linearized 
strains.  The  power  residue  method  for  generating  pseudo- 
random numbers  is  employed  in  the  technique  presented  for 
constructing  the  random  function,  whose  statistical  proper- 
ties correspond  closely  to  those  of  pressure  signals 
measured  in  the  noise  field  of  a turbulent  subsonic  air  jet. 

The  nonlinear  equation  of  motion  for  the  beam  is 
written  in  finite  difference  form  and  programmed  together 
with  the  random  load  simulation  in  FORTRAN  (709)  for 
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i 

solution  on  a digital  computer.  The  forcing  function  used 
represents  a time-random  concentrated  load  applied  trans- 
versely at  midspan.  Boundary  conditions  used  include 
rotationally  fixed  and  pinned  ends. 

Numerical  solutions  are  obtained,  from  which  sta- 
tistical measures  of  response  at  midspan  and  at  the  quarter 
points  of  the  beam  are  computed.  Response  autocorrelation 
functions,  power  spectra,  and  probability  density  functions 
are  presented  graphically  for  a range  of  load  mean-square 
values.  The  variation  of  mean-square  load  with  mean  square 
response  is  also  shown  and  compared  with  the  corresponding 
linear  relation. 

Although  the  technique  used  in  the  present  study 
was  applied  to  the  problem  of  a beam  carrying  a time-random 
concentrated  load,  it  can  in  principle  be  applied  to  more 
complicated  structural  components  and  random  load  configu- 
rations . 


xiii 


CHAPTER  I 


INTRODUCTION 

During  the  last  decade,  increasing  attention  has 
been  given  to  the  problem  of  predicting  the  response  of 
structures  to  random  excitation  (1,2).  Practical  moti- 
vation for  such  research  has  arisen  primarily  in  the  field 
of  aircraft,  missile,  and  space  vehicle  development. 

In  the  powered  phase  of  missile  or  launch- vehicle 
flight,  both  the  missile  structure  and  payload  are  subjected 
to  random  excitation  of  high  intensity.  The  response  of 
the  total  system  to  such  loads  can  be  a determining  factor 
in  the  success  or  failure  of  a mission.  Vibration- 
sensitive  equipment,  when  carried  aloft,  can  fail  to 

t 

perform  properly  if  the  vibrational  environment  is  too 
severe;  over  a sufficiently  long  operating  life,  random 
excitations  can  contribute  to  the  fatigue  failure  of  a 
structural  component. 

The  earliest  analyses  of  the  response  of  structures 
to  random  loads  were  concerned  with  linear  structures,  i.e. 
those  having  linear  differential  equations  of  motion  (1). 
Physical  structures,  however,  are  not  always  so  obliging 
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as  to  respond  to  excitation  in  a manner  consistent  with 
linear  theory. 

In  particular,  continuous  elastic  structural  ele- 
ments, e.g.  beams,  plates,  and  shells,  can  exhibit  large 
(finite)  deflections,  in  which  case  the  governing  equations 
of  motion  may  be  nonlinear. 

An  approximate  technique,  known  as  the  method  of 
"equivalent  linearization”  has  been  developed  and  applied 
to  several  mechanical  problems  (3).  The  method  is 
inherently  limited  to  the  class  of  problems  having  small 
nonlinearities. 

An  exact  technique  (4-)  is  available  for  determining 
an  integral  representation  of  a nonlinear  structural 
response  when  the  excitation  can  be  represented  as  a Markov 
process.  This  approach  has  been  applied  to  beams  and 
plates  with  nonlinearities  arising  from  finite,  moderately 
large  slopes  of  the  middle  surface  (5).  The  restriction 
that  an  applied  load  must  be  characterized  as  a Markov 
process  prohibits  the  use  of  this  technique  for  random 
loads  that  are  continuous  functions  of  spacial  coordinates 
or  time.  Successive  values  of  a Markov  process  are  sta- 
tistically uncorrelated,  implying  a physically  unrealizable 
power  spectrum  (white  noise)  and  an  infinite  mean  square. 
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Because  analytical  approaches  for  predicting  the 
response  of  nonlinear  structures  to  correlated  random  loads 
are  totally  lacking,  methods  of  numerical  analysis  were 
employed  in  the  present  study  to  simulate  both  a structural 
element  and  a continuous  (correlated)  random  forcing 
function. 

Numerical  techniques  for  solving  partial  differential 
equations,  both  linear  and  nonlinear  in  form,  have  been 
available  for  many  years  (6).  The  application  of  such 

i 

techniques  has  of  late  been  rendered  practicable  by  the 
existence  of  high  speed  digital  computers.  But  before 
such  techniques  can  be  employed  in  a study  of  random 
vibrations,  it  is  necessary  that  a sequence  of  random 
numbers,  representing  a desired  forcing  function,  be  made 
available.  The  present  study  was  therefore  initially 
concerned  with  the  development  of  an  efficient  means  of 
simulating  physically  realizable  random  loads  in  a form 
directly  usable  in  a digital  computer.  (That  method  of 
simulation  is  described  in  Chapter  II.)  The  structural 
model  chosen  for  investigation  was  that  of  a beam  whose 
resistance  to  deformation  arises  from  bending  and  stretch- 
ing. The  response  of  such  a beam  to  transverse  loading, 
for  moderately  large  slopes  of  the  middle  surface,  is 
linear  when  the  ends  are  unrestrained  against  longitudinal 
translation  and  nonlinear  when  both  ends  sire  restrained 


against  longitudinal  translation.  Both  cases  were  con- 
sidered, with  the  forcing  function  in  the  form  of  a 
concentrated,  time-random  load  at  midspan,  and  boundary 
conditions  corresponding  to  rotationally  unrestrained 
ends  and  rotationally  fixed  ends. 


CHAPTER  II 


SIMULATION  OF  A RANDOM  FORCING  FUNCTION 

Physical  loads  associated  with  continuous  structures 
are  most  generally  functions  of  time  and  spatial  coordi- 
nates. For  simplicity,  this  chapter  is  concerned  with  the 
simulation  of  a continuous,  Gaussian,  random  function,  F, 
of  one  variable,  t.  The  method  employed  can,  by  direct 
extension,  be  applied  to  simulate  continuous,  Gaussian, 
random  functions  of  a finite  number  of  variables. 

Before  proceeding  to  a description  of  the  random 
load  simulation,  those  statistical  and  related  concepts 
basic  to  the  analysis  will  be  considered  briefly. 

2.1.  Random  numbers  and  Pseudo-Random  Numbers 

Under  its  classical  definition,  a set  of  random 
numbers  is  distinguished  not  so  much  by  its  properties  as 
by  its  historical  origin  (7).  Random  numbers  are  numbers 
assigned  to  the  elements  of  a sample  space  which  in  turn 
are  defined  as  the  possible  outcomes  of  some  chance  event 
or  experiment  (8).  For  example,  a gambler’s  die  has  six 
faces  distinguished  by  dots  corresponding  to  the  numbers 
one  through  six.  If  the  die  is  true  and  "fairly”  tossed, 
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each  face  has  an  equal  probability  of  being  on  top  when  the 
die  comes  to  rest.  (The  probability  that  a particular 
outcome  will  occur  in  a given  experiment  is  the  number  of 
results  corresponding  to  that  outcome  divided  by  the  total 
number  of  possible  results.)  Thus,  for  a perfect  die 
given  a fair  toss,  the  number  of  dots  on  the  top  face,  when 
the  die  comes  to  rest,  is  by  definition  a random  number. 
Random  numbers  then  can  be  generated  by  gambling  devices 
(to  the  extent  that  such  devices  are  unbiased)  or  by 
physical  processes  such  as  radioactive  decay.  It’is  thus 
possible  to  generate  and  store,  in  printed  tables  or  on 
punched  cards,  large  sequences  of  random  numbers. 

As  a practical  matter,  however,  it  is  often  incon- 
venient to  use  such  tables  or  cards  when  high  speed  digital 
computers  are  employed.  A less  time  consuming,  and  there- 
fore more  efficient,  procedure  is  for  the  computer  to 
generate,  deterministically,  a set  of  numbers  which  are 
operationally  random,  i.e.  which  are  indistinguishable 
from  random  numbers  under  statistical  tests  of  randomness. 
Such  numbers  have  been  given  the  name  "pseudo-random 
numbers'  (9)»  Hereafter  we  shall  no  longer  be  concerned 
with  random  numbers,  but  only  with  pseudo— random  numbers. 

To  avoid  a cumbersome  phrase,  pseudo-random  numbers  will 


henceforth  be  referred  to  simply  as  random  numbers. 


7 


Many  methods  have  been  devised  for  generating  such 
random  numbers.  One  of  the  most  successful,  called  the 
power  residue  method  (9,10),  was  employed  in  the  present 
study. 

The  following  is,  in  brief,  the  power  residue  method. 
Let  SQ  be  any  odd  integer.  Let 

(2.1)  S/  = K S(._, 

If  b is  the  word  size  of  a binary  computer,  K is  chosen 
such  that  K ^ 2*^  , Calculation  of  KSQ  produces  a 
product  ft  bits  long.  The  multiplication,  performed  with 
double  precision,  permits  the  b least  significant  digits 
to  be  retained  as  the  next  random  integer,  while  the  b 
higher  bits  are  discarded.  The  sequence  of  random 
integers  generated  by  continuation  of  this  procedure  are 
converted  to  a set  of  random  numbers  (called  the  "basic 
set")  which  are  distributed  uniformly  over  the  unit 
interval  bounded  from  below  by  zero,  by  floating  the 
random  integers  and  establishing  a decimal  point  to  the 
immediate  left  of  each  integer. 

The  sequence  of  numbers  thus  obtained  will  eventually 
repeat  itself,  i.e.  the  sequence  is  periodic.  However, 
the  35  hit  word  length  of  the  IBM709  computer  used  in  this 
study  permits  a sequence  with  a period  of  "over  8.5 
billion  numbers"  to  be  generated  (10).  Thus,  the 
periodicity  of  those  numbers  was  ignored. 
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2.2  Statistical  Background. 

Let  V(t)  be  a function  generated  by  a random  or 

"chance"  process.  The  value  of  the  function  at  any 

t «=  t is  thus  a random  variable,  -oo<v<  +oo.  Associated 
o 

with  v is  a function  f(v),  which  for  v continuous  is 
defined  such  that 

(2.2)  £(v)  ^ Oj 

(2.3)  f{  (v)dv  = i , 

— oo 

rb 

(2A)  J f(v)d\r  = p{  3 - V ^ b } J 
a 

where  a and  b are  any  two  values  of  v , with  b > a , 
and  p{a^v^bj'  is  the  probability  that  v lies  within 
the  closed  interval  bounded  by  a and  b . The  function 
f(v)  is  called  the  probability  density  function  of  v . 

The  mean  of  v , denoted  by  fjv  is  defined  as 

(2.5)  = J v f rv)  dv  . 

• oo 

The  mean  is  a measure  of  the  central  tendency  of  v . 

The  r a central  moment  of  v is  defined  as 

(2.6)  jJv  = j C^-^)r  f(v)dv. 

- OO 


The  averages  indicated  in  definitions  (2.5)  and 

(2.6)  are  designated  ensemble  averages.  Of  particular 
statistical  interest  is  the  second  central  moment  (r«2)» 
commonly  called  the  variance,  which  is  a measure  of  dis- 
persion about  the  mean.  The  square  root  of  the  variance 
is  designated  standard  deviation.  When  - 0 in 
expression  (2.6),  the  variance  is  called  the  mean  square. 

The  function  V(t)  may  also  be  averaged  over  t . 

If  V(t)  is  continuous  and  defined  for  - oo  < t < +oo  , 
then 

(2.7)  = lim  - f VCt)  Jt  , 

I Tl  — — iT 

(2.8)  Vrv  = lim  ^T[Wt) 

' I T|  — > °°  J 

and,  in  particular, 

(2.9)  = ul  = lim  ~r  f [V(t)  - pv]  dt  j 

I t|  *■  °°  J_r 

where  a bar  indicates  averaging  over  the  range  of  t . 

A function,  V(t),  is  said  to  be  stationary  if 
averages  over  t , with  /T/  finite  and  sufficiently 
large,  are  independent  of  the  location  of  the  origin  t * 0 
A stationary  random  function  may  also  be  ergodic. 

A random  function,  V(t),  for  which  averages  over  the  range 
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of  t equal  the  corresponding  ensemble  averages,  is  said 
to  be  ergodic.  The  significance  of  ergodicity  is  that  any 
member  of  the  ensemble  of  an  ergodic  process  is  statistical- 
ly typical  of  the  ensemble.  (It  is  of  value  to  note  that 
every  random  process  that  is  stationary  and  has  zero  mean, 
Gaussian  probability  density  function,  and  a continuous 
autocorrelation  function  is  also  ergodic.) 

The  autocorrelation  function  of  V(t)  is  defined 
as 

(2.10)  Rv(rj  = lim  ±f  V(t)  l/ct  +t;  dt  . 

J Tl  — *■  oo  J 

By  letting  z - t -T  it  can  be  shown  that  R v (~T)  = 


i.e.  Rvfr)  is  an  even  function.  It  should  be  noted  that 

(2.11)  Hv(o)  = lira  ±f ; 

Jt|  —*oo  J _ 


that  is,  Ry  Co)  is  the  mean  square  of  V(t).  In  the 
absence  of  periodic  or  constant  components  of  V(t), 

(2.12)  lim  RV(~T)  = o. 

171  — 

The  power  spectral  density  function  or  power 
spectrum  is  defined  as  the  Fourier  transform  of  the  auto- 
correlation function: 

1 

See  Reference  2,  p.  10. 
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(2.13)  5 v ( co)  — __ _ 


/ 


The  power  spectrum  Sv (oo)  associates  a component  of  the 
mean  square  of  V(t)  with  every  real  value  of  the  fre- 
quency go  . 

The  relation  between  an  ergodic  random  function 
and  its  autocorrelation  function  is  established  as 
follows.  (A  random  function,  W(t),  defined  over  a finite 
range  of  t will  be  considered  because,  in  practice,  only 
finite  records  of  random  functions  are  encountered.)  Let 
W(t)  represent  a sample  of  an  £rgodic  random  process  such 
that  W(t)  » 0 for  t ^ 0 and  for  t ^ T , The  Fourier 
transform  of  W(t)  is 


Sv(cu,T)  is  a random  variable  that  is  a function  of  T 
and  takes  on  different  values  for  different  members  of 
its  ensemble.  In  addition, 


O 


O 

The  power  spectrum  of  W(t)  is  then 


(2.16)  lim  ^[SffCLOjT)]  - Sw(tu) 


T-» 


See  Reference  2,  p.  12 
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where  £ indicates  expected  value  and  the  limit  indicated 
is  the  same  for  all  members  of  the  ensemble.  The  autocor- 
relation function,  a3  has  been  noted,  is  the  Fourier 
transform  of  S^(gj,T).  Thus,  if  two  ergodic  random 
functions  of  finite  duration  are  typical  samples  of  the 
same  ensemble,  their  autocorrelation  functions  will  be  the 
same  within  a small  variation.  Or,  for  such  ergodic 
functions  having  the  same  probability  density  function, 
the  more  nearly  coincident  the  autocorrelation  functions, 
the  higher  the  percent  statistical  confidence  with  which 
it  may  be  said  that  the  ergodic  functions  belong  to  the 
same  ensemble. 

2«3  Gaussian  Random  Humber  Generation 

In  the  present  study,  the  power  residue  method 
discussed  in  section  (2.1),  above,  was  programmed  for  use 
on  an  IBM  709  digital  computer.  For  reference,  that 
program  is  included  in  the  Appendix,  where  it  appears  in 
the  Subroutine  for  Uncorrelated  Standard  Normal  Variable. 

A statistical  test  to  provide  a measure  of  the 
uniformity  of  the  distribution  of  the  basic  set  was  per- 
formed. Ideally,  the  elements  of  the  basic  set,  denoted 
by  , should  be  distributed  uniformly  over  a unit 
interval,  indicating  an  equal  probability  of  occurrence 
for  each  of  the  B^.  Accordingly,  1124  elements  of  the 


13 


basic  set  were  generated  and  classified  into  20  equal  sub- 
intervals  of  the  unit  interval.  The  expected  number  of 
elements  in  each  subinterval  was,  accordingly  1124  -f  20  - 
56.2.  Observed  frequencies  of  occurrence,  shown  in  Table 
1,  are  clustered  about  the  expected  value.  A quantitative 
measure  of  the  uniformity  of  the  distribution  is  given  by 
the  linear  correlation  coefficient,  R , of  the  statistical 
frequencies  shown  in  column  2 of  Table  1.  An  R » 0 


corresponds  to  a zero  slope  for  the  linear  regression  line 
and  hence  is  associated  with  a uniform  distribution.  (R 

o 

is  defined  as 


where  - the  number  of  the  classification  interval,  ot^  ■ 
the  frequency  of  occurrence  in  the  Kib  interval,  Z is  the 
mean  of  the  classification  intervals,  and  N,  in  the  present 
case,  is  20.)  The  calculated  value  of  RQ  for  the  was 

found  to  be  -0.000985,  indicating  an  apparently  negligible 
slope  for  the  regression  line.  This  implies  that  the  Bi 
are  uniformly  distributed  to  a very  good  approximation. 

The  elements,  B^,  of  the  basic  set  were  employed 
to  construct  a set  of  random  numbers  having  a Gaussian 


o 


N 


(2.1?)  R0  = J<£. 

N 


Z fz  -zKk 

K =1 


DISTRIBUTION  OP  1124  ELEMENTS  OP  THE  "BASIC  SET 
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probability  density  function,  as  follows.  The  probability 
density  function  of  the  being 


(2.i8)  fra;  - 


'l  , * ^ b ^ l 

o } B <o  ; B > i , 


the  population  mean,  fJ 3,  and  variance,  , are 

(2.19)  IJB  = f BJB  = ~ , 

0 

r ? 2. 

(2.20)  -y  (8  - <JB  = ~ . 


Let  Bj  denote  elements  of  the  set  of  sample  means  ob- 
tained from  samples  of  size  n taken  from  set  B.  That  is, 


J + n 


(2.21)  Bj  = Bj  j j = 

i*  J 

By  the  central  limit  theorem  of  statistics,^  the  b\  have 
a probability  density  function  that  approaches  a Gaussian 


distribution  as  n 


Further,  the  mean,  | Jg,  of  the 


— 1 2 P 

population  B is  jj  B ^ and  the  variance,  ^ , is  /n 

• • A compact  expression  for  the  previous  statement 

is  that  the  Bj  are  distributed  as  N(^,^^),  indicating 

a normal  (Gaussian)  distribution  with  mean  and 


3 


See,  for  example,  Reference  8,  p.  107. 
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variance  . It  follows  from  the  above  that  if  a finite, 
but  sufficiently  large,  n is  chosen,  the  distribution  of 


the  Bj  will  approximate 


The  distribution  of  the  B.  differs  from  that  of  an 

J 

ideally  Gaussian  random  variable  in  that  the  range  of  the 

former  is  (0,1)  while  the  range  of  the  latter  is  (-oo,  +00), 

It  is  therefore  necessary  to  choose  n such  that  the 
2 

variance  cr-  is  small  enough  to  render  the  probability 

of  occurrence  of  a B^  ^ 1 or  a Bj  < 0 negligible. 

The  number  of  standard  deviations  Ng  between  pg  and 

the  extreme  values  B.  - 0 , B.  - 1 is 

0 J 


(2.22)  Ns  = 1---^  = = _J_ 

B Og  zldl 

/ /2n 


= vnr  . 


For  n - 25,  N-  * 8.66.  The  probability  of  achieving  a 
value  of  B^  more  than  8.66  standard  deviations  from  jj  ^ 
is  less  than  0.00006.  Therefore  the  possibility  of 
achieving  B^  > 1 or  B^  < 0 was  considered  negligible; 
the  value  n - 25  was  employed  in  the  present  investigation 
and  was  found  to  be  satisfactory. 

The  above  method  of  generating  values  of  a Gaussian 
random  variable  by  sampling  from  a uniformly  distributed 
random  variable  was  programmed  for  solution  on  a digital 
computer  and  appears,  together  with  the  power  residue 
method,  in  the  Appendix,  where  it  appears  in  the  Subroutine 
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for  Uncorrelated  Standard  Normal  Variable.  Each  entry 
into  the  subroutine  of  ENTRY  STDNR  causes  25  values  of  B. 
(called  LSTN  in  the  program)  to  be  generated  by  the  power 
residue  method.  From  those  values,  a Bj  is  computed  by 
equation  (2.21),  above,  and  made  available  within  the 
main  program  under  the  symbol  XX. 

A^ocorr elated  Random  Function  Simulation 

From  the  set  of  numbers  B,  a continuous,  piecewise 
linear  function,  Q(t),  is  constructed  as  follows. 

Pairs  of  Gaussian,  random  numbers, 

(2.23)  - Bj  and 

(2.24)  Tt  -yT  Bj+1  , i . 1,2, ...q  , j . 1,3,5,...  r, 

are  defined  in  terms  of  the  elements  of  set  B.  Here  U 
, __  ’ r T 

and  c^T  are  the  preassigned  mean  and  standard  deviation 

of  the  T±  , which  are  interpreted  as  random  increments  of 

a continuous  variable  t . That  is,  t m . £ T where 
. . 1-1  1 ’ 

ie  a particular  value  of  t . Corresponding  to  t m 

is  a value  of  Q(t)  denoted  by  and  determined  from 

expression  (2.23).  By  linearly  interpolating  between  suc- 
cessive pairs  , t1 ; . . . iQ^t^,  a continuous,  piecewise 
linear  function,  Q(t),  0 « t * T,  is  obtained.  (Because 
a successful  simulation  was  obtained  with  linear  interpo- 
lation, no  more  complicated  interpolating  functions  were 
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investigated.)  A continuous,  piecewise  linear,  random 

function,  F(t),  having  a predetermined  mean,  pv,  and 

■ — 2 

variance,  f ±s  obtained  from  Q(t)  as  follows.  Let 

QCt)  +<*  j o ^ t ^ T 

(2.25)  F(t)  = < 

[ O J ? <o  j <>r  j 

or 

(2.26)  Fj  * J3Qt  +«,  J = l,  ~jn. 

Then, 

(2.27)  pF  =-hL  Fr 

A ^ “*  / 

(2.28)  ^‘=]j-£[F/-pF]S  = y^  . 

/ 

Thus , ^ and  a - - J3  fJ  ^ , whence 

(2.29)  - “rf"  [ Q(  t)  + pp  - jJQ  ] 

In  the  present  investigation,  pF  is  chosen  to  be  zero, 
whence 

(2.50)  J(t)  - [ Q(  O ~ p„J  . 

The  computer  program  for  constructing  F(t)  bj  the  method 
outlined  above  appears  in  the  Appendix  as  Segment  1 of  the 
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main  program.  There,  the  interpolated  values,  F^,  , of 
F(t)  are  available  for  equal  increments  of  time,  S t, 
under  the  FORTRAN  symbol  FSMALL  (I). 

The  percent  statistical  confidence  with  which  it 
may  be  said  that  F(t)  has  a Gaussian  probability  density 
function  was  measured  directly  using  a chi-square  test. 

The  quantity  X2  is  defined  as 


where  0^  is  an  "observed"  value  of  a random  variable, 
ei  is  the  "expected"  value  of  that  variable  and  m is 
the  number  of  intervals  into  which  the  range  of  the  random 
variable  is,  for  convenience,  divided.  A value  of  X2  =>  0 
indicates  that  0^  = , i • l,...,m  and  thus  that  the 

observed  values  may  be  said  to  have  a Gaussian  distribu- 
tion with  100  per  cent  statistical  confidence. 

To  test  the  statistical  hypothesis  that  F(t)  , 
represented  by  the  discrete  interpolated  values  F^  , has 

a Gaussian  distribution,  the  population  mean,  and 

p * 

variance,  p , of  F(t)  must  be  known;  as  predetermined 

measures  of  F(t),  both  pp  and  <rp2  are  known  quanti- 
ties . 

The  range  of  t , i.e.  JJ  + 8.66  was  divided 
into  23  intervals  to  provide  a chi-square  test  with  23 


statistical  degrees  of  freedom.  In  any  interval,  i, 
bounded  by  values  Fa  and  (a  > b),  the  expected 

number  of  values  of  F^  is 

fF< 

(2.32)  et  = nj  0(F)  JF, 

n 

where  n is  the  number  of  values  of  F^  upon  which  J? 
is  to  be  based,  and  0(F)  is  the  Gaussian  probability 
density  function, 

(2-55)  0(F^ 

Integrals  of  the  form  given  by  expression  (2.32)  were 
approximated  by  summations 

20 

(2.34)  e } « , +F  “ziCFi-F^j 

Jml 

and  programmed  for  solution  on  a digital  computer.  That 
routine  appears  as  Segment  2 of  the  program  given  in  the 
Appendix.  As  a check  on  the  accuracy  of  the  numerical 
integration,  a test  quantity  corresponding  to 


(2.35) 
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was  computed  and  given  the  symbol  UNITY,  (if  the  range  of 
Integration  for  integral  (2.35)  is  extended  to  (—,  + oo)> 
the  value  of  that  integral  is  unity  in  accordance  with 
expression  (2.3)  of  the  definition  for  a probability 
density  function.)  The  calculated  value  of  UNITY  was 

1.000000  indicating  an  overall  accuracy  of  at  least  seven 
significant  figures. 

The  value  of  X2,  as  calculated  from  expression 
(2.31)  and  (2.32)  based  upon  5000  interpolated  values  j , 
was  found  to  be  3.496,  which  for  23  degrees  of  freedom 
implies  that  F(t)  may  be  said  to  belong  to  a Saussian 
distribution  with  better  than  98  per  cent  statistical 
confidence,^ 

For  a stationary  F(t),  the  normalized  autocorrela- 
tion function,  By,  is  a function  of  y only.  (The  auto- 
correlation function,  expression  (2.10),  is  normalized  by 
dividing  by  the  mean  square,  5,(0).)  That  is,  finite 
samples  F(t),  0-t  ^ T,  will,  for  T sufficiently  large, 
yield  the  same  Rp  to  a good  approximation. 

Sp  was  calculated  from  the  F,  using  numerical 
integration,  based  upon  the  series  approximation 

_ Xj  Fj/Fj+K 

(2.36)  , K‘T/tt, 


See  Reference  8,  Table  III,  p.  313. 


as  programmed  in  Segment  3 of  the  Appendix.  The  range,  n 
of  is  a measure  of  the  averaging  time,  T * n fit  1 
being  an  even  function  of  tT  , calculations  were  made  only 
for  non-negative  values  of  y . Functions  Ip  , for  F(t) 
are  shown  in  Figure  1,  for  averaging  times  ranging  from 
5C0  fit  to  4500  fit.  The  curves  of  Figure  1 indicate  that 
for  averaging  times  T ^ 1000  fit,  Rp  behaves  in  a nearly 
stationary  manner,  the  variations  from  the  curve  for 
being  ^ 0.03  in  absolute  value. 

The  particular  value  of  fit  to  be  chosen  must  be 
such  that  solutions  of  the  equation  of  motion  discussed 
below  are  numerically  stable  and  convergent.  In  addition 
the  time  scale  thus  established  should  be  consistent  with 
a realistic  Rp  . Both  criteria  can  be  satisfied  since 
8 c • X-j/n  where  fit  is  chosen  to  equal  At  (the  finite 
increment  of  t for  which  solutions  mentioned  above  are 
numerically  stable  and  convergent),  Tj  , is  the  value  of  T 
at  the  first  zero  crossing  of  an  empirically  determined 
autocorrelation  function,  lp,  and  n is  the  value  of  the 
abscissa  at  the  first  zero  crossing  in  the  curve  for 
T/  fit  - 1000  in  Figure  1.  The  value  of  n is  in  turn 
dependent  upon  ^ , the  preassigned  average  number  of  F^  ®s 
in  an  interval  J^r  . (Although  no  rational  functional  rela- 
tion between  n and  § is  known  to  the  author,  it  was 
found  in  practice  that  their  magnitudes 


are  very  nearly 
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equal  for  the  values  of  n (or  $ ) used  in  the  present 
study.)  The  empirical  autocorrelation  function,  1 (T  ) 
used  is  that  recorded  by  T^ubert,  Kizner,  and  Nash  (11,12). 
The  works  cited  describe  statistically  the  pressure  signals 
obtained  in  the  noise  field  of  a turbulent  subsonic  air 
jet. 

Figure  2 provides  a comparison  between  R (T)  and 
the  empirical  autocorrelation  function,  I ( T ) , mentioned 
above.  The  curve  Rp(T  ) is  there  plotted  with  n . 78  , 
from  Figure  1 and  hence  St  ~ 0.6  x 10"5  seconds;  the 
effect  is  to  cause  Rg  to  have  the  same  first  zero  cross- 
ing point  as  R-g  . The  following  should  be  noted:  (1) 

Rj,  and  Ip  are  indistinguishable  for  /'T / - 0.5  milli- 
seconds, (2)  the  absolute  areal  difference  between  the 
curvas  of  Rp  and  R^  is  0.0784  milliseconds,  i.e.  10.4 
per  cent  of  the  gross  area  under  the  empirical  curve,  and 
(3)  the  form  of  Rj,  conforms  approximately  to  that  of  Rg  . 
In  the  frequency  domain,  the  power  spectrum  of  F(t)  has 
an  apparent  upper  bound  at  approximately  1300  cycles  per 
second  as  does  that  for  Fp(t).  The  mean  frequency  for 
the  bandwidth  of  F(t)  is  600  c.p.s.  as  compared  with  a 
mean  frequency  of  550  c.p.s.  for  the  same  bandwidth  of  ' 


8*0 
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The  function  F(t),  generated  by  the  method  intro- 
duced in  this  chapter,  was  shown  above  to  be  (with  better 
than  98%  statistical  confidence)  a Gaussian  random  function. 
For  averaging  times  T 2 1000  ft,  F(t)  was  shown  to  be 
stationary.  The  mean  of  F(t)  was  predetermined  to  be 
zero.  And  F(t)  being  continuous  by  definition,  its 
autocorrelation  function  is  also  continuous.  Thus,  having 
satisfied  the  criteria  noted  in  section  (2.2),  F(t)  has 
been  shown  to  be  ergodic.  (Results  obtained  below  by 
using  F(t)  as  a forcing  function  therefore  apply  to  an 
entire  class  of  random  loads  — that  class  having  the  same 
statistical  properties  as  the  sample  F( t ).) 

The  empirical  random  function,  Fg(t) , of  which 
is  the  autocorrelation  function,  is  also,  to  a good 
approximation,  Gaussian  and  stationary,  has  a zero  mean 
and  continuous  autocorrelation  function,  and  hence  is  also 
ergodic  (12).  As  is  shown  in  Figure  2,  Rp  closely 
approximates  Kg  , with  St  - 0.6  x 10^  seConds.  F and 
F£  both  being  ergodic,  and  their  autocorrelation  functions 
differing  only  by  small  variations,  the  ensembles  of  which 

F and  Fg,  respectively,  are  members  have  similar  sta- 
tistical properties. 


CHAPTER  III 


MATHEMATICAL  MODEL  OF 
THE  EFFECTS  OF  FINITE 

OF 


THE  BERNOULLI -EULER  BEAM  INCLUDING 

lateral  deflection  and  stretching 

THE  MIDDLE  SURFACE 


In  linear  beam  theory,  a Bernoulli -Euler  beam  is 
one  in  which  the  resistance  to  applied  loads  is  assumed 
to  arise  from  a bending  couple  proportionate  to  the 
curvature  of  the  bent  beam.1  However,  if  the  beam  ends 
are  restrained  from  translational  displacements,  the 
consequent  equation  of  motion  is  exactly  linear  only  for 
infinitesimal  deflections.  The  nonlinearity  of  the 
equation,  for  finite  deflections,  arises  from  the  presence 
of  a term  associated  with  the  stretching  of  the  middle 
surface  of  the  beam.  The  assumptions  underlying  such  a 
nonlinear  3ernoulli-Euler  beam  model  are  presented  in  the 
following  section.  The  beam  model  itself  is  a special 
case  of  the  well  known  Von  Karman  plate  equations  (13). 

Ur,  Equation  of  Motion  Under  Transverse  Load 

The  beam  under  consideration  is  referred  to 
Lagrangian  cartesian  coordinates  x,y,z.  The  middle 
surface  of  the  beam  is  initially  coincident  with  the  x-z 

1 — 

See  Reference  13,  pp.  2-3. 
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plane,  the  beam  being  oriented  along  the  positive  x axis 

Motion  is  restricted  such  that  displacements  are  parallel 
to  the  x-y  plane. 

Neglecting  longitudinal  inertia  and  damping,  the 
differential  equations  of  motion  become 


(5.1) 

(5.2) 


Jx 


= o J 


Lfldydx  - OTi  dy  £ w ( cr$  -f  l J*  )/  b tv  , 

by  J bx  { *X  A Jx 

lj"alx)c/y  +l2*/clxJy  + 3 Jy  <Jx  1 w _ 

***  ^ 'by  * ix 


ft  Ijft  o/xc/  y j 

^ bt*  J 


where  ~ represents  deflection  in  the  y direction  and^o 
is  mass  per  unit  area  parallel  to  the  x-y  plane. 

Integrating  equation  (3.1)  across  the  depth,  h,  of 
the  beam  yields 


(3.3) 

Letting 


A 

2 


/ azJy  + °y  J = 


2 


O. 


h 

2 


A 

•2 


-h 

2 


y 


and  assuming  zero  shear  stresses  at  the  surfaces  of  the 
beam,  equation  (3.3)  becomes 

(5.4)  i_A/,  = 
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Thus,  at  any  instant  of  time,  Nx  is  constant  throughout 
the  length  of  the  beam. 

Equation  (3.2)  can  be  rewritten  as 

4-  t °*y  -i. 

TIT 


(3.5)  -h  q-  i«v 


y*U/ 


0 w _ /) 

} •V'  y 


where  a term  of  order  greater  than  two  has  been  neglected. 

Integrating  equation  (3.5)  across  the  depth  of  the  beam 
yields 


(3.6)  op  I 


+ tS ^ * 4v  <?*  + ^f*-A/x  + 


:fc 

J 


if  / -/ 


Aa! 


i X { bx 


'di2 


where 

Assuming 


2 ^vo/, 


*J  °7  * 


<7 


L 


- c 


L 

y = ~2 


where  <j  represents  a lateral  surface  load,  and  recalling 
that  the  term  in  braces  is  zero,  equation  (3.6)  becomes 

(3.7)  q + Ai*  ^ £)  = oAAv 

I %XZ  dx  * -X 
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Multiplying  equation  (3.1)  by  y and  integrating 
over  the  depth  of  the  beam  yields 


(3.8)  >_./**  + [ ]_f 

(3.9)  LM..-0  , 


- Q. 


= o 


where  M - 
x 


1 

2 


'L 

2 


Noting  from  equation  (3.9)  that 
(3.10)  £ Q - 


ax 


a*' 


and  substituting  equation  (3.10)  into  equation  (3.7)  yields 

(3.n)  <?  + N*  -f.  ilMx  = />A  a2w 

' ax2  ^ a 

Quantities  Mx  and  Nx  can  be  expressed  as  functions  of 
derivatives  of  deflection  w , Assuming 


it  follows  that 


(3.13) 


/■ 


<•« 

c/5  - L 


E 


where 


s 


is  measured  along  the  longitudinal  axis  of  the 
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deformed  beam,  a is  the  cross-sectional  area  of  the 
beam,  E is  Young's  modulus,  and  h is  the  length  of  the 
deformed  longitudinal  axis.  Noting  that 


(5.14)  c/s  = X/x*  tJw1  = yT  -A  (iff)1  Jx 


it  follows  that 


(5.15)  Js  = r ! . J/ijV  _ , , 

L + * (w  s(h<j 1 J* , 

(5.16)  J>S  a £ 


Substituting  equation  (3.16)  into  equation  (3.13)  yields 

a(.  -(¥/(»' j„. 

<>  ^ 

Let  q consist  of  a forcing  function  $"(x,t) 
and  a linear  damping  term  . aen)  substituting 

equation  (3.17)  and  the  moment-curvature  relation, 

(3.18)  M*  = -El  £"  a 

into  equation  (3.11),  the  equation  of  motion  is  obtained 
in  terms  of  deflection  w : 
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(3.19)  El 


ZL  >X*  L zxj 


bX+ 
d 2 tv 


- SA\g  K 


bx 


'iW 

TE 


+ J-  (*>t)  + 


whore  y3  is  now  considered  to  be  mass  per  unit  volume. 

The  difficulty  of  obtaining  analytical  solutions 
for  the  nonlinear  integro-differential  equation  (3.19) 
is  no  obstacle  to  the  numerical  solution  of  that  equation, 
assuming  j"(x,t)  to  be  known.  The  following  finite  dif- 


ference analogs  are  employed: 


(5.20) 

V'V'  (V,  t ) 

'■‘V.y 

(5.21) 

) IV 

ax 

2 (AX) 

(3.22) 

^tV 

H4>/,y  -2-wc\j  + wt- 

_/.»y 

(3.23) 

ax2 

.4 

d w 

a x * 

4 fflV6y 

(3.24) 

o/  VvY  . 

Woj+» 

" rnx;  4 

o It 

2 (n't; 

where  i and  j correspond,  respectively  to  the  spatial 
coordinate,  x,  and  time,  t.  Symbols  Ax  and  At  denote 
finite  increments  of  x and  t , each  of  which  is 


sufficiently  small  to  assure  convergence  and  stability  of 
the  numerical  solution. 
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Substituting  relations  (3.20)  through  (3.23)  into 
equation  (3.19) ♦ and  letting 

(3.25)  Y,-.j  = , 


results  in  the  following  equations  of  motion: 


(3.26) 

(3.2?) 


=.  y cj  j 

ol~t 


and 


p/  y<  = — — — -4 

o j-t 

VS/  t + 6 j 


L 


IV 


<+zJ 


- 4 f IV, 


J 


~ W c-  z jJ  J 


JL  y 

">*  '"J 


+ 


f*  y + slIm3  w<»  + 

N 

^ <-'>j ) i 


where  stations  i * 2 and  i » N define  the  locations  of 
the  beam  ends  x = 0 and  x * L , respectively. 

Substituting  relation  (3.24)  and  a similar  one  for 
into  equations  (3.26)  and  (3.27)  yields 

(3.28)  w<tj+i  - 2 (a±)  Yi, j •+  wvv^y-/  j 
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i(4t)EJ  r 
(AX)+  L 


£ C 

+ Lp(<xJ*  1 IV'*V 


Note  that  a separable  solution  has  not  been  assumed. 

Together  with  boundary  and  initial  conditions,  equations 

(3.28)  and  (3.29)  constitute  an  algorithm  for  determining 
tV(x,t)  numerically. 

1^.2  Boundary  and  Initial  Condition* 

Physical  boundary  conditions  may  impose  varying 
degrees  of  restraint  against  rotation  and  translation  of 
the  beam  ends.  The  two  limiting  cases  of  rotational 
restraint  pinned  and  fixed  ends,  are  considered  in  the 
present  study.  Similarly,  the  cases  of  no  translational 
restraint  and  complete  translational  restraint  at  the  ends 


are  considered 
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A pinned  constraint,  by  definition,  offers  no  re- 
sistance to  rotation.  The  corresponding  analytical 
boundary  conditions  are 

(3»30)  WCcLj't)  ~ o } 

(3.31)  MxC*,'t)  = oJ 

where  a denotes  the  value  of  x at  a beam  end.  In  view 
of  equation  (3.18),  condition  (3.31)  may  also  be  expressed 
as 

(3.32)  = o. 

lxz 

In  finite  difference  form,  letting  m denote  an 
end  station  of  the  beam  (m  =*  2,N),  pinned  boundary  con- 
ditions (3.30)  and  (3.32)  become 

(3.33)  =0, 

(3.34)  . 

A fixed  end  constraint,  by  definition,  permits  no 
rotational  motion  of  a beam  end.  The  corresponding 
analytical  boundary  conditions  are,  therefore, 

(3.35)  wfij-O  = o, 

b W Cdj  i)  = 

ax 


(3.36) 
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In  finite  difference  form,  fixed  boundary  conditions 
(3.35)  and  (3.36)  become 

(3.37)  VV/h^y  »Oj 

(3*38)  +/} j ~ W m~  1 3j  . 

The  initial  conditions  to  be  imposed  are  quite  arbi- 
trarily selected.  The  effect  of  initial  conditions  on 
the  response  of  a damped  beam  will  subside  after  a suf- 
ficiently great  lapse  of  time.  In  the  present  study,  sta- 
tistical averages  were  taken  over  increasing  averaging 
times  until  a stationary  response  was  obtained.  The 
simplest  initial  conditions,  those  of  static  equilibrium, 
were  employed: 

(3.39)  va /(XjO)  = o3 

(3.40)  (X'°)  =r  o 

it 

In  finite  difference  form,  initial  conditions  (3.39)  and 
(3.^0)  become 

(3.41)  =o 

(3.42)  tV/jt  j 


where  j - 1 corresponds  to  t « 0 . 
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The  equations  of  motion  (3.26)  and  (3.27),  together 
with  the  boundary  and  initial  conditions  discussed  in  this 
section,  were  programmed  for  solution  on  a digital  computer. 
The  program  appears  as  Segment  4 of  the  main  program  in  the 
Appendix.  The  forcing  function  employed  was  a time-random 
concentrated  load,  P(t),  located  at  midspan: 

(3.43)  = Pit)  S(X- L/d), 


where  N - 6, 10,1-4-, 18,  ••  • . 

The  range  of  N was  chosen  for  convenience  in  de- 
termining the  response  at  midspan  and  at  the  quarter 
points  of  the  beam. 


and  P(t)  - (Ax)P(t)  . 


In  finite  difference  terms, 


CHAPTER  IV 

DEFLECTION  RESPONSE  OF  BEAMS  TO  DETERMINISTIC 
AND  CORRELATED  RANDOM  FORCING  FUNCTIONS 


In  solving  a diff erential  equation  by  numerical 
methods,  criteria  must  be  established  for  determining  the 
stability  and  convergence  of  solutions.  For  all  but  the 
simplest  equations,  no  general  criteria  are  known. ^ In 
their  absence,  the  most  common  technique  used  is  to  re- 
peatedly reduce  the  magnitude  of  the  finite  increments 
until  no  appreciable  change  is  observed  in  the  response. 
That  technique  was  employed  in  the  present  study. 


4.1  Deflection  Response  for  a Concentrated 
Magnitude  Suddenly  Applied  at  Midspan 


Load  of  Constant 


To  determine  the  range  of  values,  for  At  and  Ax, 
in  which  stable  and  convergent  numerical  solutions  of 


equation  (3.19)  might  be  expected  for  a random  forcing 
function,  a preliminary  study  was  made  of  a beam  excited 
by  a constant  magnitude  concentrated  load  applied  suddenly 
at  midspan.  The  end  boundaries  were  assumed  to  be  pinned 
and  restrained  against  translational  motion.  Damping  was 


T 


See  Reference  14,  p.  489. 
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taken  to  be  zero.  The  values  of  beam  parameters  used  are 
listed  in  Table  2.  They  correspond  to  a rectangular  steel 
beam  of  unit  width,  with  a depth-to-length  ratio  of  1/20. 

Many  different  values  for  s and  x were  employed  in 
successive  integrations  of  the  equation  of  motion.  It  was 
found  that  numerical  instability  of  the  solution  occurred 
for  values  of  At  :>  0.02  milliseconds.  The  response  was 
found  to  converge  well  for  Ax<L/16  « 2.5  and  At  ^0.006 
milliseconds . 

Deflections  at  the  midspan  (x=L/2)  and  quarter 
points  (X-L/4.3LA)  of  the  beam  are  shown  in  Figures  3 
and  4 where  P(t)  was  taken  as 

f i^O 

(4.1)  F ]Ci)  = 1 

[ 0 , ±<0 

and 


(4.2)  /jet; 


12-jSOO*6 j t 

o j t < O 


respectively.  The  magnitudes  of  and  P2  were  chosen 

to  obtain  deflections  well  within  the  linear  range  of 

response  for  Px  and  well  within  the  nonlinear  range  of 
response  for  P2« 

The  curves  of  Figures  3 and  4 suggest  a predominantly 
first  mode  response;  w(L/2,t),  w(LA,t),  and  w(3L/4,t)  are 
m phase  and  the  deflection  vs.  time  curves  are  sinusoidal 


TABLE  2 
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numerical  values  op  beam  parameters 


Symbol 

Numerical  Value 
(Units:  in. -lb.- 
— , sec,  system') 

L 

40 

w 

1 

h 

2 

E 

30  x 106 

/> 

.732  x 10“^ 

l/2)  w(x  = L/4,3L/4) 


<}• 


OJ  r-\  O 


03 

Tj 


d 


o 

o 

03 

a 


r-t 

a 

d 

•H 

43 

•H 


Bsqoti-p  cp  -[00*0  x uof409-[j©(i 
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Pig.  3. — Deflection  response  for  a concentrated  load  of  25  lb.  applied 
suddenly  at  midspan 


L/2)  w(x  - L/4,3L/4) 
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Pig.  4a — Nonlinear  deflection  response 
applied  suddenly  at  midspan 


4-3 

in  appearance.  For  purposes  of  comparison,  the  linear 
static  deflections  that  would  obtain  at  midspan  under 
slowly  applied  concentrated  loads  of  the  same  magnitude  as 
the  step  functions,  P^  and  P2,  were  computed.  Those  calcu- 
lated deflections  were  based  upon  an  equivalent  linear 
spring  constant,  Kv,  defined  as 


(4-. 3)  ke 


P _ 49  El 
IA/  L3 


The  amplitudes  of  the  curves  in  Figure  3 indicate  that  the 
corresponding  deflections  can  be  anticipated  to  lie  well 
within  the  linear  range  of  response.  (The  maximum  de- 
flection-to-depth-of-beam  vartio  is  0.0017.)  The  double 
amplitude  of  w(L/2,t)  should  therefore  be  exactly  twice 
the  static  deflection  mentioned  above.  This  expectation  is 

in  fact  borne  outs  the  double  amplitude  of  w(L/2,t)  is 

_2 

0.3334-  x 10  inches  while  the  corresponding  static  de- 

-2 

flection  is  0.1666...X  10  inches.  Thus,  with  an 
accuracy  of  (0.3334-  x 10“^  - 2 x 0.1666...X  10"*^)  (0.333 

...x  10  ) - 0.021  per  cent,  the  double  amplitude  of 

w(L/2,t)  equals  twice  the  corresponding  static  deflection. 

In  Figure  4-,  the  effects  of  restraining  the  beam 
ends  against  translation  are  observable.  The  frequency  of 
the  response  has  increased,  as  compared  with  the  very 
nearly  linear  case  shown  in  Figure  3.  The  linear  static 


deflection  at  midspan  is,  for  Figure  4,  approximately  0.83 
inches.  Thus,  the  double  amplitude  relative  to  the 
corresponding  linear  static  deflection  has  been  reduced  as 
compared  with  the  linear  case  illustrated  in  Figure  3. 
Stretching  of  the  middle  surface,  caused  by  restraining 
the  ends  against  translation,  is  sufficient  to  account 
for  both  the  increase  in  frequency  and  relative  decrease 
in  amplitude. 

Statistical  Response  .to  a Time-Random  Concentrated 
Load  Applied  at  Midsoan  ~~  ' 

The  equation  of  motion  (3.19)  in  finite  difference 
form  (equations  3.28  and  3.29)  was  solved  using  the  random 
load  F(t)  as  a forcing  function.  Numerical  solutions  were 
carried  out  for  several  combinations  of  boundary  conditions 
and  preassigned  values  of  o'2 . Values  of  beam  parameters 
used  are  those  listed  in  Table  2.  In  addition,  the 
damping  coefficient,  K,  was  taken  to  be  0.5. 

Samples  of  deflection  responses,  measured  at  midspan 
and  at  the  quarter  points  of  the  beam,  are  shown  for  each 
of  the  cases  studied  in  Figures  5 through  16.  Response 
means  and  variances  were  computed  under  Segment  5 of  the 
main  program  given  in  the  Appendix.  In  no  case  was  a 
response  mean  greater  than  55  microinches  obtained,  indi- 
cating that  all  response  means  were  found  to  be  very  nearly 


V2  x « L/4,3L/4 
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0.4  r x - L/2  x « 1/4,31/4 
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L/2  x * L/4,3L/4 
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L/2  x = L/4,3L/4 
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L/2  x « L/4,3L/4 
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55 


N 

X 


• ••*••••  • 

ooooooooo 

I I I I 


CO 

4 


vo 

KN 


4 

OJ 


OJ 


o 


w 

d 

o 

o 

© 

CQ 

•H 


a 

d 


0) 

a 

•H 

EH 


BsqotiT  rtf  uoj^oexjsa 


•d 

© 


X 

•rH 

m 


t) 

d 


© 


s 

o 

P. 

m 

© 

P 

d 

o 

•H 

-P 

O 

© 


© 

T) 

s 

© 

d 


0 

*? 

1 

LT\ 
i — f 


bO 

•H 


L/2  x = L/4,3L/4 


56 


57 


equal  to  zero.  Response  standard  deviations  are  compared 
with  the  corresponding  standard  deviations  of  the  excita- 
tion in  Figures  17  and  18.  Statistical  frequency  function 
approximations  of  the  response  probability  density  functions 
were  obtained  under  the  same  segment  of  the  main  program 
and  are  shown  for  the  linear  cases  and  those  nonlinear 
cases  differing  appreciably  from  the  normal  distribution 
in  Figures  19  through  24. 

Response  autocorrelation  functions  and  power  spectra 
were  computed  under  Segment  6 of  the  main  program  (see 
Appendix)  and  are  shown,  for  the  cases  studied,  in  Figures 
25  through  48.  The  autocorrelation  functions  are  shown  in 
normalized  form  to  facilitate  comparison.  Corresponding 
power  spectra  shown  are  based  upon  the  un-normalized 
autocorrelation  functions. 

Linear  cases  were  established  by  introducing  the 
constraint  Nx  » 0,  which  corresponds  to  beam  ends  unres- 
trained against  longitudinal  translation. 

Numerical  solutions  of  equations  5.28  and  3.29  were 
apparently  convergent  for  Ax  £ 2.5  inches  (L/16)  and  At  £ 0.6 
x 10  seconds.  Various  averaging  times  were  used  in  com- 
puting the  response  statistics.  Stationary  output  was 
obtained  in  all  cases  for  averaging  times  T^/At  ^ 8000. 

The  most  salient  feature  of  the  deflection  vs.  time 
curves  of  Figures  5 through  16  is  that  in  all  cases,  / 
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Fig.  18. — Excitation  and  response  standard  deviations  (ends  fixed) 
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Pig.  20. — Nonlinear-response  probability 
ends  pinned) 
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cr(x  - L/2)  - 1.038;  cr(x  > L/4,3L/4)  - 0.733 
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64 


•=t 


ends  fixed) 


<r(x  - L/2)  - 0.453;  <r Cx  - L/4,3L/4)  « 0.248 


65 


ends  fixed 


66 


noi^BX8aJ:oooq.nB  pezfrenuoii 


Fig.  25. — Nonlinear-response  autocorrelation  functions  (o'  - 1,00c)1 
ends  pinned)  ^ 
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Fig.  27. — Nonlinear-respouse  autocorrelation  functions  (cr  * 10,000; 
ends  pinned)  P 
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fixed) 
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Fig*  35* — Nonlinear-response  autocorrelation  functions  (o'  * 20,000' 
ends  fixed)  p 
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Fig.  37. — Nonlinear-response  power  spectra  (cr  - 1,000  , ends  pinned) 
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Fig.  59. — Nonlinear-response  power  spectra  (01  ■ 10,000  , ends  pinned) 
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Fig.  43. — Nonlinear-response  power  spectra  (cr  « 1,000  , ends  fixed) 
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Fig.  44. — Nonlinear-response  power  spectra  (or  « 5,000%  ends  fixed) 
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Fig.  47. — Nonlinear-response  power  spectra  (cr  - 20,000^,  ends  fixed) 
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w(LA,t),  w(L/2,t),  and  w(3L/4,t)  are  "in  phase"  in  the 
sense  that  w(L/4,t)  and  w(3LA,t)  correspond  approximately 
to  the  deflection  v(L/2,t)  at  a reduced  ordinate  scale. 
This  feature  of  the  curves  is  consistent  with  a predomi- 
nantly first  mode  response.  It  is  of  interest  to  note 
that  Herbert  (5),  in  considering  the  response  of  a similar 
beam  to  uncorrelated  (white  noise)  excitation,  found  that 
the  first  mode  response  was  a good  approximation  to  the 
total  response,  although  higher  modes  contributed  to  the 
first  mode  itself  (the  modes  being  coupled  in  the  non- 
linear case).  Figures  1 7 and  18  indicate  that  departures 
from  linear  response  (ends  unrestrained  against  longitudi- 
nal translation)  are  more  significant  for  the  more  flexible 

pinned-end  beam  than  for  the  stiffer  rotationally  fixed- 
end  beam. 

The  nonlinear-response  probability  density  functions 
(Figures  19  through  24)  exhibit  a characteristic  increase 
in  the  probability  density  in  the  neighborhood  of  the  mean 
deflection  (w.O)  and  a corresponding  decrease  in  the 
probability  density  for  more  extreme  values,  both  effects 
being  considered  relative  to  a Gaussian  distribution. 

This  effect  is  consistent  with  the  tendency  of  the  stretched 
middle  surface  to  restore  the  beam  to  equilibrium  from  any 
deformed  position.  The  "distortion"  of  the  probability 
density  functions  from  the  Gaussian  form  is  not  appreciable 
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for  cr^  * l,OOO^f  5,000#  in.  the  pinned  end  case  and  for 
crp  « 1,000#,  5,000#,  10,000#,  and  15,000#  in  the  rota- 
tionally  fixed-end  case. 

The  response  autocorrelation  functions  (Figures  25 
through  56)  exhibit  the  characteristic  form  that  is  des- 
cribed qualitatively  as  a variation  of  a damped  cosine 
curve.  Physical  interpretations  of  these  curves  can  be 
more  easily  visualized  by  considering  them  in  the  frequency 
domain,  that  is,  by  considering  their  corresponding  power 
spectra  (Figures  37  through  48). 

The  most  salient  feature  of  the  power  spectra  is, 
in  all  cases,  the  presence  of  a single  prominent  peak. 

For  increasing  cr^y  two  flanking  peaks  of  considerably  less 
magnitude  emerge.  The  frequency  responses  shown  indicate 
that  the  beam  has,  in  effect,  acted  as  a band  pass  filter 
with  regard  to  the  input  power  spectrum.  The  similarity 
in  form  of  the  response  power  spectra  for  x = L/2  and 
x « L/4 , 3L/4  and  the  prominent  peak  in  each  case  is  con- 
sistent with  the  predominantly  first  mode  response  affected 
by  higher  modes  mentioned  by  Herbert  (5). 


CHAPTER  V 


CONCLUSIONS 

The  principal  conclusions  to  be  derived  from  the 
present  study  are  as  follows: 

1)  Using  the  method  described  in  Chapter  II,  above,  it 
is  possible  to  generate  a continuous,  ergodic,  random 
function  whose  statistical  properties  closely  approximate 
those  of  pressure  signals  obtained  experimentally  from 
the  noise  field  of  a turbulent  subsonic  air  Jet. 

2)  Numerically  stable  and  convergent  solutions  for  the 
finite  difference  equivalents  of  the  equation  of  motion 
of  a beam,  whose  resistance  to  deformation  is  due  to 
bending  and  stretching,  can  be  obtained  with  numerical 
values  of  Ax  and  At  large  enough  to  render  the  use  of 
existing  digital  computers  practicable. 

3)  Solutions  for  a time-random  concentrated  load  applied 
transversely  at  midspan  contain  strong  indications  that 
the  total  response  is  largely  associated  with  the  first 
mode,  whose  frequency  varies  with  the  contribution  of 
higher  modes  in  the  nonlinear  case. 
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4)  For  a given  mean-square  load,  mean-square  response 
was  found  to  depart  more  from  the  linear  response  for  the 
relatively  flexible  pinned-end  beam  than  for  the  stiff er 
fixed-end  beam. 

Finally  it  should  be  noted  that  although  the 
technique  presented  in  this  investigation  was  applied  to 
a beam  carrying  a concentrated  load,  the  method  is  in 
principle  applicable  to  more  complicated  structural 
components  and  random  load  configurations. 


APPENDIX 
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DIGITAL  COMPUTES  PROGRAM:  'SIMULATION; 

RANDOM  VIBRATIONS  OF  BEAMS' 


1 READ  INPUT  TAPE  5 , 301 , FS , SMDT , NOD , LAP , NUMB , NT , NCELL 

301  FORMAT  (2F7. 0,215,313) 

READ  INPUT  TAPE  5 , 302 , EL , VI , D , RHO , CEE , E , AA 

302  FORMAT  (3F4.1 ,4F9.1) 

WRITE  OUTPUT  TAPE  6, 303, FS , SMDT, NOD, LAP, NUMB, 

1 NT, NCELL, EL, VI, D, RHO, CEE, E,AA 

303  FORMAT  (1H1/2B15.5,5IS///7E15.5) 

DIMENSION  FSMALL  (10000) ,BETA(23) ,GAMMA(23) , 

1 GG(400)  ,C0R(400)  ,PA(40)  ,P(40)  ,RA(40)  ,R(40) ,QUAG(40)  , 
1 EA(40) ,Y(40) ,V(40)  ,<*A(2000)  ,03(2000)  ,QC(2000) , 

1 FF (2000 ) , PRZ ( 50 ) , PHI ( 1 500 ) 

SEGMENT  1.  GENERATION  OF  RANDOM  FORCING  FUNCTION 


XX  - 

T - 

TT  - 
SFSM 


0.0 

0.0 

0.0 

- 0.0 


10 


11 


12 

13 

102 


14 


SFSMS  - 0.0 
I - 0 
H - SMDT 
CALL  RESET 
DTMU  - AA*SMDT 
SIGDT  - 0. 35*  DTMU 
FPEAKA  - 0.0 
CALL  STDNR(XX) 

FPEAKB  - XX 
NCHK  - 0 
CALL  STDNR(XX) 

DTPEAK  - DTMU  + SIGDT* XX 
I F ( DTPEAK -SMDT ) 1 2 , 1 4 , 1 4 
NCHK  - NCHK  + 1 
IF  (NCHK  - 50)11,11,13 
WRITE  OUTPUT  TAPE  6.102 

FORMAT  (1H0/1H0/30X20HEXECUTI0N  TERMINATED/// 

1 15X50H50  CONSECUTIVE  VALUES  OF  DTPEAK  ARE  LESS 
1 THAN  SMDT) 

GO  TO  1 

SLOPE  » ( FPEAKB-FPEAKA ) /DTPEAK 

I - I + 1 

FSMALL(I)  * FPEAKA  + SLOPE* (T-TT) 

TT  - TT  + DTPEAK 
GO  TO  18 
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15  T a T + SMDT 

IF  (TT-T)16, 16,17 

16  FPEAKA  - FPEAKB 
GO  TO  10 

17  JA.  - I 

I - I + 1 

FSMALL(I)  = FSMALL(JA)  + SLOPE* SMDT 

18  SFSM  - SFSM  + FSMALL(I) 

3FSMS  - SFSMS  + FSMALL(I)**2 
NFSM  - NFSM  + 1 
IF(NFSM-NOD)15,19,19 

19  TOTFSM  = I 

FSMBAR  « SFSM/TOTFSM 

FSMBSq  - FSMBAR* *2 

SSUBFS  - SFSMS/TOTFSM-FSMBSq 

22  FSTDEV  . SqRTF( SSUBFS) 

Tq  - FS/FSTDE V 

DO  200  I - 1,N0D 

200  FSMALL(I)  a Tq* (FSMALL( I) -FSMBAR) 

END  OF  SEGMENT  1. 

SEGMENT  2.  CHI  SqUARE  TEST  TO  DETERMINE  THE 
PER  CENT  STATISTICAL  CONFIDENCE  WITH  WHICH  THE 
RANDOM  FORCING  FUNCTION  MAY  BE  SAID  TO  BE  GAUSSIAN. 
DDD  - -8.66025*FS 
WW  - 17.3205*FS/FL0ATF(NCELL) 

WWW  = WW/10. 

EEE  = 2. *FS*  *2 
GGG  - DDD 
BOT  - DDD 
TOP  - BOT  + WW 
EPS  - 0. 

D025K  a 1 ,NCELL 
ALPHA  - 0. 

23  ALPHA  a ALPHA  + WWW/EXPF(GGG* *2/EEE) 

GGG  * GGG  + WWW 

IF  (GGG-TQP) 23 ,24,24 

24  BETA(K)  a FLOATF ( NOD )* ALPHA *0. 3989^2 3/FS 

EPS  - EPS  + BETA(K) 

BOT  - BOT  + WW 

25  TOP  . BOT  + WW 

UNITY  = EPS/FLOATF(NOD) 

D036J  - 1 ,NOD 

BOT  - DDD 

TOP  a BOT  + WW 

D035K  - l.NCELL 

IF  (FSMALL(J)-B0T)34,52,31 

31  IF  (T0P-FSMALL(J))34,32,33 

32  GAMMA(K)  a GAMMA(K)  +0.5 
GO  TO  36 


33  SAMMA(K)  - GAMMA(K)  +1.0 
GO  TO  36 

34  BOT  - BOT  + WW 

35  TOP  - BOT  + W 

36  CONTINUE 
PAH  = 0 . 

D038K  - 1,NCELL 

38  PAR  - PAR  + (GAMMA(K)-BETA(K) )*  *2/BBTA(K) 
WRITE  OUTPUT  TAPE  6, 106, PAR, UNITY 
106  FORMAT ( 1H1/1HO/25X , 5HCH1SQ20X , 5HUNITY/ 

1 20X , FI 2 . 6 , 1 3X , FI 2 . 6// / ) 

END  OF  SEGMENT  2. 


39 


40 


45 

50 


107 

1 


108 

60 


SEGMENT  3.  AUTOCORRELATION  FUNCTION 

D039I  « 1,NT 

GG(I)  s 0,0 

JBB  - NOD-NT 

LL  - 1 

SQ  - 0.0 

D060L  - LAP, JBB, LAP 

D040J  - LL,L 

SQ  - SQ  + FSMALL(J)**2 

D050-I  - 1 ,NT 

D045J  - LL,L 

K « I + J 

GG(I)  * GG(I)  + FSMALL(J)*FSMALL(K) 

COR(I)  - GG(I)/SQ 
LL  - L + 1 

WRITE  OUTPUT  TAPE  6,107,L 

FORMAT  (1H1/IH0/20X,4HL  - I4//14X.1HI23X. 

GHCOR( I )////) 


WRITE  OUTPUT  TAPE  6,108, (I ,COR(I) ,1*1 ,NT) 
FORMAT( 10X , 15 , 10X , E20 . 7 ) 

CONTINUE 

END  OF  SEGMENT  3. 


SEGMENT  4.  NUMERICAL  INTEGRATION  OF  EQUATION 
OF  MOTION. 

MM  = NUMB-2 
MB  - 1 
MO  * MM/4 
NALPHA  = MO  + 2 
NBETA  - MM/2  + 2 
NGAMMA  » NUMB-MO 
KEGA  =.  NUMB  + 1 
KEGB  = NUMB  - 1 
GO  » MM 
DX  ■ EL/GO 
A - WI*D 

BAR  - (WI*D**3)/12. 


POE  - l./(RHO*A) 

CPOE  - CEE*PQE 

VA  - E*BAR*POE/ ( DX*  *4) 

VB  - VA*A*DX*0.125/(BAR*EL) 

HB  . 2. *H 
D0300I  - 1,KEGA 
PA(I)  = 0.0 
P(I)  - 0.0 
RA(I)  » 0.0 
300  R(I ) - 0.0 

D0500J  - 1 .NOD 
QUAG(NBETA)  » FSMALL(J) 

ZIG  - 0.0 
D0350I  - 2, NUMB 
K - I + 1 
L - I - 1 

350  ZIG  * ZIG  + (P('K)-P(L) ) * *2 
D0400I  = 3,KEG3 
1*1  + 1 
M * I + 2 
L » I - 1 
N - I - 2 

EA(I)  - -CP0E*PA(I)-VA*(P(M)-4.*(PCK)+-P(L)) 

1 +6.*P(I)+P(N))+VB*(P(K)-2.*P(I)+P(L))*ZIG 
1 +POE*QUAG(I) 

Y(I)  - RA(I)+HB*3A(I) 

400  W(I)  «*  R(I)+HB*PA(I) 

W(KEGA)  = W(KEGB) 

W(l)  * W(3) 

Y(KEGA)  » Y(KEGB) 

Y(l)  =>  Y(3) 

MA  - MA  + 1 
IF  (MA-5)415, 410,410 
410  MA  = 0 

QA(MB)  - W(NALPHA) 
qB(MB)  « W(NBETA) 

QC(MB)  - W(NGAMMA) 

MB  = MB  + 1 
415  D0420I  » l.KEGA 

RA( I ) - PA(I) 

R(D  - P(I) 

PA(I)  « Y(I) 

420  P(I)  - V(I) 

500  CONTINUE 

WRITE  OUTPUT  TAPE  6,501 

501  FORMAT( 1H1/1H0/20X , 9HW( NALPHA) 10X , 8HW(NBETA) 
1 12X,3HW(NGAMMA)///) 

WRITE  OUTPUT  TAPE  6,502, (^A(J) ,QB(J) ,QC(J) , 

1 J - 20, MB, 20) 

502  F0RMAT(10X,3E20.5) 

END  OF  SEGMENT  4. 
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SEGMENT  5.  MEAN , VARIANCE , AND  PROBABILITY 
DENSITY  FUNCTION  OF  RESPONSE. 

D0510I  = 1 ,MB 
RT  =»  RT  + QA(I) 

RS  = RS  + QB(I) 

RTZ  = RTZ  + QA(I)**2 
510  RSZ  = RSZ  + QB(I)**2 
FMB  = MB 
QAMEAN  * RT/FMB 
QBMEAN  - RS/FMB 
QAUANC  = RTZ /FMB  - QAMEAN * * 2 
QBUANC  - RSZ/FMB  - QBMEAN* *2 
WRITE  OUTPUT  TAPE  6 , 512 , QAMEAN , QBMEAN , 

1 QAUANC. QBUANC 

512  FORMAT(lHl/lHO/33X,9HQAMEAN  « E10.5,10X, 

1 9HQBMEAN  = ElO. 5,//55X,9HQAUANC  - E10.5,10X, 
1 9HQBUANC  = ElO. 5////) 

D0550J  - 1,MB 

550  FF( J)  - QA(J) 

FMN  = QAMEAN 
FU  - QAUANC 

WV  =,  0.12*SQRTF(FU) 

DDD  - FMN-3.  *3QRTF(FU) 

GO  TO  600 
552  D055H  - 1,50 

551  PRZ(I)  - 0.0 
DO555J  « 1 .MB 

355  FF(J)  « QB(J) 

FMN  - QBMEAN 

FU  . QBUANC 

WW  - 0. 12*SQRTF(FU) 

DDD  =■  FMN-3. *SQRTF(FU) 

600  D0650J  - 1,MB 

BOT  » DDD 
TOP  - BOT  + WW 
D0635K  -1.50 
IF(FF(J)-B0T)634, 632,631 

631  IF(TOP-??(J))634,632,633 

632  PRZ(K)  * PRZ(K)  +0.5 
GO  TO  650 

635  PRZ(X)  = PRZ(K)  +1.0 
GO  TO  650 

634  BOT  - BOT  + WW 

635  TOP  « BOT  + WW 

650  CONTINUE 

DO  62 OK  * 1,50 

620  PRZ(K)  - PRZ(K)/(FMB*WW) 

IR  - IR  + 1 

IF( IR-2 ) 660 , 670 , 670 


660  WRITE  OUTPUT  TAPE  6,700 

700  F0RMAT(40X,41HPR0BABILITY  DENSITY  FUNCTION 
1 OF  V(N ALPHA )////) 

WRITE  OUTPUT  TAPE  6,701 , (I ,PRZ(I) ,1  - 1,50) 

701  F0RMAT(40X,I2,E20.5) 

GO  TO  552 

670  WRITE  OUTPUT  TAPE  6,702 

702  FORMAT ( IH1/1H0/40X , 40HPR0BABILITY  DENSITY 
1 FUNCTION  OF  W(NBETA)///) 

WRITE  OUTPUT  TAPE  6,701 , (I ,PRZ(I) ,1=1,50) 

END  OF  SEGMENT  5. 

SEGMENT  6.  RESPONSE  AUTOCORRELATION  FUNCTIONS 
AND  POWER  SPECTRA. 

DG750I  * l.MB 
750  FF(I)  - QA(I) 

GO  TO  800 

755  DO  7601  - 1,MB 
760  FF(I)  - Q3(I) 

800  DO  8391  » 1 ,NT 
COR(I)  = 0.0 

839  GG(I)  - 0.0 
JBB  - MB-NT 
LL  - 1 

SQ  - 0.0 
MBS  - MB/10 
MBA  * 9*MBS 
D0860L  * MBA, MB, MBS 
D0840J  * LL,L 

840  SQ  > SQ  + FF( J) * *2 
D0850I  = 1 ,NT 
D0845J  » LL,L 

K - I + J 

845  GG(I)  » GG(I)  + FF( J) *FF(K) 

850  COR(I)  « GG(I)/SQ 
LL  - L + 1 

WRITE  OUTPUT  TAPE  6,900,L 
900  F0RMAT(1H1/1H0/20X,11H0UTPUTL  * I4//14X, 

1 1HI23X,6HC0R( I )////) 

WRITE  OUTPUT  TAPE  6,908, (I,COR(I) ,1=10, NT, 10) 
908  F0RMAT(10X,I5,10X,E20.7) 

D 08411  - 1 ,NT 

841  COR(I)  - GG( I ) 

DTAU  - 5. *SMDT 
DOUB  = 2.* DTAU 
DARG  - 6. 28318* DTAU 
D0855I  “ 50,1500,50 
SLAM  =0.0 

D0854K  - 1,NT 
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l 


t 


854 

855 

910 

911 

860 


870 

909 


c 


STDNR 


1) 


XR 

RESET 


GEN 

CONS 

LSTN 

SET 

SUM 


SLAM  . SLAM  + COR(K) *C0SF(FL0ATF(I*K) *DARG) 

PH1(I)  * D0UB*(0.5+SLAM) 

WRITE  OUTPUT  TAPE  6,910 

FORMAT ( ///l OX , 39HSPECTRALDENSITYCORRESPONDING 
1 TO  ABOVE/I OX, 24H AUTOCORRELATION  FUNCTION//) 

WRITE  OUTPUT  TAPE  6,911 
FORMAT ( 1 4X , 1 H J2  3X , 6HPHI ( J ) // ) 

WRITE  OUTPUT  TAPE  6, 908, (I, PHI(I) ,1-50,1500,50) 
CONTINUE 
NCOB  - NCOB  + 1 
IF(NCOB-2)755 . 870 , 870 
WRITE  OUTPUT  TAPE  6,909 
FORMAT( //10X , 41HEND  OF  RESPONSE 
1 AUTOCORRELATION  FUNCTIONS) 

END  OF  SEGMENT  6. 

CALL  EXIT 

END 

FAP 

SUBROUTINE  FOR  UNCORRELATED  STANDARD  NORMAL  VARIABLE 

ENTRY  STDNR 

ENTRY  RESET 

SXA  XR,1 

STZ  SUM 

AXT  25,1 

LDQ,  GEN 

MPY  LSTN 

STQ  LSTN 

CLA  LSTN 

ARS  8 

ORA  CONS 

FSB  CONS 

ACL  -1B8 

FAD  SUM 

STO  SUM 

TIX  1),1,1 

CLA  SUM 

FDH  - 25. 

XCA 

FSB  CONS 
XCA 

FMP  . 17.5205 
STO*  1,4 
AXT  **,1 
TRA  2,4 
CLA  SET 
STO  LSTN 
TRA  1,4 

DEC  1220703125 
DEC  0.5 
BCI  1,C00001 
BCI  1,000001 

END 
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